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Abstract. We study several dynamical properties of a recently proposed 
implementation of the quantum transverse-held Ising chain in the framework of circuit 
QED. Particular emphasis is placed on the effects of disorder on the nonequilibrium 
behavior of the system. We show that small amounts of fabrication-induced disorder 
in the system parameters do not jeopardize the observation of previously-predicted 
phenomena. Based on a numerical extraction of the mean free path of the system, 
we also provide a simple quantitative estimate for certain disorder effects on the 
nonequilibrium dynamics of the circuit QED quantum simulator. We discuss the 
transition from weak to strong disorder, characterized by the onset of Anderson 
localization of the system's wave functions, and the qualitatively different dynamics it 
leads to. 
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1. Introduction 

Circuit quantum electrodynamics (QED) systems consist of superconducting artificial 
atoms coupled to the electromagnetic field in a microwave resonator [Tj. Such 
systems have been successfully used for implementations of elementary quantum optical 
Hamiltonians [2] [3] and basic quantum information processing [U IH [6] [7] . The rapid 
technological development in the field of circuit QED will soon facilitate experiments 
with highly coherent multi-atom, multi-resonator circuit QED architectures. This makes 
circuit QED a promising platform for observing interesting multi-atom quantum optical 
effects [El [9] [TO] and even for simulating genuinely interacting quantum many-body 
systems from solid state physics HU [121 [131 HH HS1 HS1 HZl HH HH [20]. 

In [20] . we have proposed and analyzed a circuit QED design that implements 
the quantum transverse-field Ising chain (TFIC) coupled to a microwave resonator for 
readout. The TFIC is an elementary example of an integrable quantum many-body 
system. Despite its simplicity, it still exhibits interesting features, e.g. a quantum phase 
transition (QPT), and therefore serves as a model example system in the theory of 
quantum criticality [21] and nonequilibrium thermodynamics [22] . Our circuit QED 
quantum simulator can be used to study quench dynamics, the propagation of localized 
excitations, and other nonequilibrium phenomena in the TFIC, based on a design that 
could easily be extended to break the integrability of the system. While in [20] we have 
focussed on an idealized implementation of the TFIC with perfectly uniform parameters, 
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the main purpose of the present article is to investigate the effects of disorder in the 
system parameters on the dynamical behavior of our quantum simulator. 

The study of disorder effects on quantum simulators is relevant for two reasons. 
First, on the more practical level, any real experimental system will come with a degree 
of unwanted disorder (especially in condensed matter settings). In the case of circuit 
QED systems, inhomogeneities of the system parameters are caused by fabrication 
issues as well as by static noise fields (e.g. produced by defects). It is important to 
verify that the basic behavior of a quantum simulator survives the amounts of disorder 
which are present in realistic systems, or even to estimate the amount of disorder that 
can be tolerated. Second, on the more fundamental level, simulating quantum many- 
body systems with built-in (potentially tunable) disorder is interesting in its own right. 
Many physical phenomena, from free propagation of wave packets to quench dynamics 
to (quantum) phase transitions can be affected in significant ways by disorder, and this 
leads to phenomena such as Anderson localization or disorder-induced phases. 

To prepare our study, we briefly review the system (section 2.1), discuss sources 
of disorder and how disorder scales with the tunable system parameters (section 2.2), 
and explain the mathematical approach to and some properties of the quantum Ising 
chain (section 2.3). We start our main discussion by considering the time-dependent 
correlations of the order parameter of the chain, where the finite-size effects and the long- 
time behaviour will be analyzed in the absence of disorder (section 3.1). Based on this, 
we will move on to the spectrum of the resonator coupled to the quantum Ising chain in 
our system, which is closely related to the aforementioned time-dependent correlations. 
To that end, we employ a very useful approximation which we have introduced in [20] 
and which will presumably become important also for future studies of quantum many- 
body systems coupled to resonators. In this approximation, the full quantum many- 
body system is replaced by a bath of harmonic oscillators with identical spectrum. 
We show here that this approximation actually works very well under appropriate 
circumstances (section [3T2 ) . We then calculate the spectrum of the resonator coupled 
to a slightly disordered Ising chain and find that the effects of disorder on the spectrum 
are small (section 3.3). The Ising chain in our circuit-QED quantum simulator can be 
driven out of equilibrium in several ways. This allows one to perform various types 
of nonequilibrium experiments, a particularly appealing application of our setup. In 
[20], we have suggested to observe the propagation of a localized excitation through the 
chain or the nonequilibrium dynamics of the system after a quantum quench. Here we 
show that the predicted phenomena are insensitive to a small amount of disorder in 
the system parameters (section 4.1 and section 4.2, respectively). However, as argued 
above, it would be highly desirable to possess also a quantitative theory of disorder 
effects. Since the nonequilibrium dynamics of the TFIC is determined by the ballistic 
propagation of quasiparticles (wave packets), we formulate and numerically verify a 
relation between the mean free path of the latter and the parameters of the system 



CONTENTS 



4 



a 



3 



N 



-resonator A / resonator B- 

■fv) „t„ m o (optional) 



artificial atom j 



Figure 1. Circuit QED implementation of the quantum transverse-field fsing chain 
(adapted from 20J ) . Charge-based artificial atoms are capacitively coupled to their 
nearest neighbors. Coupling the first (iVth) artificial atom to resonator A (B) allows 
one to use standard circuit QED techniques for initialization and read-out of the first 
(iVth) artificial atom. 



and the disorder potential. By means of this relation we are able to differentiate weak 
from strong disorder, to predict the dynamical behavior of our quantum simulator given 
a certain disorder strength, and to estimate the amount of disorder that a particular 
experiment can tolerate (section 4.1). 



2. The quantum transverse-field Ising chain in circuit QED 

2.1. Setup 

We consider a circuit QED quantum simulator of the TFIC as proposed in [20J. It 
consists of a chain of iV capacitively coupled charge-based superconducting artificial 
atoms [23], such as transmons or Cooper-pair boxes (the latter have to be biased to 
their charge degeneracy point [23J to properly simulate the TFIC [20] )• For a review on 
superconducting artificial atoms, see [23] . The first artificial atom is capacitively coupled 
to a microwave resonator (see Figure [T]). This resonator A is required for initialization 
and readout of the first artificial atom. For certain types of experiments, e.g., for 
measuring end-to-end correlators, one also needs a second resonator B, coupled to the 
iVth artificial atom. For details on the implementation and the theoretical description of 
the system, see [20J. The system (at first, only with resonator A) can be approximately 
described by the Bamiltonian 

H = u cra + g(al + a)al + Hi, (1) 

and Hi is the Bamiltonian of the TFIC, 

N N-1 

% = E^-E^ ( 2 ) 

3=1 3=1 

Bere, a^, is a Pauli matrix. That is, the artificial atoms are considered as two-level 
systems (qubits), and their two states are described as spin states. The operators a) 
and a generate and annihilate a photon of energy oj . The transition frequency Qj > 
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of the jth qubit corresponds to a local magnetic field acting on the jth spin in the 
usual interpretation of the TFIC As such, it would be transverse to the direction of 
the qubit-qubit coupling Jj. The latter can be either ferromagnetic {Jj > 0, as in the 
geometry of Figure [TJ or anti-ferromagnetic {Jj < 0, if the qubits in Figure [l] are rotated 
by 90°). While we have previously focussed in [20] on the uniform case Jj = J and 
Qj = Q for all j , we are here often interested in the case where these system parameters 
are explicitly nonuniform. This is because, on the one hand, a slight nonuniformity of 
the Qj and Jj has to be expected from imperfections of the fabrication process. On 
the other hand, one can also intentionally detune one or several qubits (by threading 
the SQUID-like loops of the qubits with different fluxes) and observe how the system's 
properties change depending on the detuning. 

2.2. Disorder and tunability of the system parameters 

Let us discuss the ffux-tunability and the undesired disorder of the system parameters 
in some more detail. We will argue that the qubit transition frequencies Qj and the 
qubit-qubit couplings Jj, when normalized to their respective mean values, may be 
assumed to be flux- in dependent. This will be relevant for our theoretical description of 
the disorder in the system. 

In reality, it should be possible to engineer the geometry of the qubits essentially 
uniform. That is, the areas of the qubits' SQUID loops, their charging energies, and the 
coupling capacitances between the qubits will only vary weakly in the chain. However, 
the (flux-tunable) total Josephson energies -Ej(<3>) of the artificial atoms should be 
experimentally harder to control since these depend exponentially on the properties 
of the Josephson junctions. For a flux-tunable (i.e. SQUID-type) artificial atom with 
two Josephson junctions |24j . 

£j (*) = (e ; + e5)»s(^)(l + ^t^(^)) 1/2 . (3) 

Here, e) is the Josephson coupling energy of one Josephson junction, <£> is the 
superconducting flux quantum, $ is the tunable external flux threading the SQUID loop, 
and d = (e] — ef)/(e] + ef). Assuming equal qubit geometries, $ can be chosen identical 
for all qubits (e.g. by using a common flux line) and only the e) can give rise to disorder. 
Even if one allows for \d\ ~ 0.1, this still means d 2 <C 1, and one can approximate the 
total Josephson energy of the jth artificial atom by Ejj{§) ps (ej - + e-}-) cos($7r/$ ) 
(as long as |$| ^ 5>o/2). Now, for Cooper-pair boxes at the charge degeneracy 
point Qj{$) ~ Ejj($), and for transmons ~ {8Ejj{$)E c ] 1/2 [211 • Thus, under 

the assumption of identical geometry, both for Cooper-pair boxes and transmons the 
transition frequencies Qj(&) of all qubits j scale with the same function a(<&) of the 
(global) flux $, Qj(®) = a{&)Qj{0). Here, a($) = cos($7r/$ ) for Cooper-pair boxes 
and a($) = [cos(<l ) 7r/ ( I>o)] 1 ^ 2 fo r transmons. This result implies that the qubit transition 
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frequencies, when normalized to their flux-dependent mean value, do not depend on $ 
and, hence, have the same statistical properties for all $. Explicitly, the mean value of 
the flj is given by £!,($) = a(§)flj(0). Thus, the mean value of the flj is flux tunable. 
However, the normalized qubit transition frequencies flj(<&)/flj(<$>) are independent of 
$, which must also be the case, for instance, for their standard deviation. This will 
become important for our numerical implementation of disorder in the flj when we 
consider changes of the external magnetic flux $. 

The qubit-qubit couplings Jj can also depend on the Ejj(&) and, thus, on the flj. 
This is the case for transmons, where approximately Jj oc (fljflj+i) 1 ! 2 oc (EjjEjj+i) 1 / 4 
[27J ED]- That is, the disorder in the flj and the Jj will not be independent for transmons. 
Moreover, flj, Jj, and their mean values flj and Jj change with the external flux $ 
approximately in the same proportion (oc [cos($7r/$ )] 1//2 ). F° r Cooper-pair boxes, on 
the other hand, the Jj depend only on charging energies and not on the Ejj(&) [20J. 
This means that the Jj are not affected by changes of the external flux. Furthermore, 
the disorder in the Jj should be less pronounced than and hardly correlated with the 
disorder in the flj. Concerning the relative strength and the correlation of the disorder in 
the Jj and the flj, we remark that also static noise fields can play a role, producing some 
disorder also in the various charging energies of the system (in particular for Cooper-pair 
boxes, which have small electrostatic capacitances). Apart from that, disorder in the Jj 
will turn out to have a much weaker effect than disorder in the flj. These deliberations 
justify to assume for simplicity that both for Cooper-pair boxes and for transmons 
disorder in the flj and Jj can be present to a comparable degree, and that disorder in 
the flj (Jj) would be uncorrelated with the disorder possibly present in the Jj (flj). We 
finally remark that many properties of the transverse-field Ising chain are determined 
by the ratio flj/ Jj, since this ratio essentially (in the limit of weak disorder) determines 
the eigenstates of the system (see below). For standard trasmons, the ratio flj/Jj is 
not straightforwardly flux-tunable. One of the experiments we suggest to do with our 
quantum simulator relies on the possibility to change the eigenfunctions of the system 



[cf. section 4.2 , which can be performed only by changing the ratio flj/Jj. All other 
possible experiments discussed in this article can be done in principle with Cooper-pair 
boxes and transmons equally well, irrespective of the Jj being flux-dependent or not [20J . 
Therefore, when plotting our results as function of a flux-tunable system parameter, we 
will assume for definiteness that our circuit QED quantum simulator of the TFIC is 
implemented with Cooper-pair boxes, and that the Jj do not change with the external 
magnetic flux. 

2.3. The transverse-field Ising chain 

The Hamiltonian ^ can be exactly diagonalized by means of a Jordan- Wigner 
transformation, which was first used in this context in [251 [26]. This transformation maps 
the spin degrees of freedom to fermionic operators Cj, cj via = cj- exp(i7r J2k=i c l c fc) 
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and yields 



N 



N 



N-l 



3 C j C 3 



i=i 



^^•[4ct +1 + C t Ci+1 +H.c. 



(4) 



Up to a constant — ^ • Qj/2, this Hamiltonian is of the form 



iV 



(5) 



Note that the conditions H = W and {c,,cj} = 1 require A = and S = —B T . 
By introducing new fermions r] k = Yl!j=\9k,j c j + hh,j c ], such Hamiltonians can be 
transformed into the diagonal form H = J^k-^kiVkVk — 1/2) + [25]. The 

components g^j and h^j of the vectors g k and and the excitation energies of H 
are determined by defining normalized vectors 4>k = gk + hk and ipk = gk — hk and by 
solving the equations 



lfc(A- J B)=A fc ^ fc , ^(A + 5)=A 



(6) 



In our case, 
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-Jn-i 



(7) 



/ 

-Ji -> Ji in A. For 



and £> is obtained by substituting Ay = Qj — y and Aj + ij — — U3 —r Uj 
uniform Qj and J}, the 0^, ^fc, and A*, can be analytically calculated from Equations 
(JgJ) (see, e.g., [20] ). For nonuniform system parameters, these quantities have to be 
determined numerically. In both cases, the Hamiltonian "Hi of the TFIC can be written 
in the form 



%-^M4-l/2), 



(8) 



and knowledge of the (f>k and ipk allows one to express spin observables in terms of the 
^-fermions, which is the basis of many of our calculations. For instance, 

a{ = (c] + Cj ){cj - ct) = E fajipk'Avl + Vk)(Vk' - 4')- ( 9 ) 



k,k' 



We collect some important facts about the TFIC. In the uniform case, 
A fc = 2jyi + £ 2 -2£cosfc. 



(10) 



Here, J = \J\ and £ = Vt/2J is the normalized transverse field. The possible values of 
are solutions of sin kN = £ sin A; (AT + 1). For A" — > oo, the uniform TFIC undergoes a 
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second order QPT at £ = ±1 from a ferromagnetic [£ G (0, 1)] or an anti-ferromagnetic 
[£ G (—1,0)] ordered phase with doubly degenerate eigenstates (one A fe — >• 0) to a 
paramagnetic disordered phase with > for all k. The QPT is signalled by the 
disappearance of long-range correlations in a x . This QPT will also occur in a nonuniform 
system (at some mean transverse field strength Qj) [2T]. However, there can be weakly 
(dis)ordered Griffith-McCoy 'phases' in the vicinity of the critical point |28" t 129 } l30 |. l3Tj . 

Finally, we introduce a convenient notation for nonuniform Qj and Jj. In this case, 
we will frequently write Qj = Qtj and Jj = j7Yj, where Tj and rj usually have mean 1, 
or, if Qj and Jj follow probability distributions, expectation value 1. We will refer to fl 
as the 'mean' qubit transition frequency, even if fl = (IX,-) is the expectation value of a 
probability distribution and the actual mean value flj is (for finite N) in general different 
from fl. We use the same convention for the qubit-qubit coupling J. Furthermore, we 
define the local and the 'mean' normalized transverse magnetic field, £j = flj/2Jj and 
£ = Q/2J. Note that in general both £ ^ £j and £ ^ (£,,■) (but for the probability 
distributions we will consider, (£ — (£?))/ (£j) < 1%)- We usually characterize %\ by the 
parameters £, J, Tj, and rj. Under the assumptions formulated in section 



2.1 



Q and 



thus £ are flux-tunable without changing the Tj in the proposed circuit QED quantum 
simulator of the TFIC 



3. Spectrum of the system 

In order to provide a guideline for the initial experimental characterization of our setup, 
we have calculated in [20] the transmission spectrum S of the resonator as a function 
of the probe frequency uj and the flux-tunable qubit transition frequency fl [see below, 



equation (28)]. To that end, we have first calculated the spectrum of the bare TFIC for 



coupling to the first qubit via a]., 

p(u) = [ dte^(^(t)^(0)), (11) 



which is the Fourier transform of the qubit autocorrelator p{t) = (o"^.(t)cr^.(0)). We have 
argued that for sufficiently large (but finite) N, qubit decay processes will render the 
measured spectrum continuous and akin to the spectrum one would obtain by taking 
the limit N — > oo in the calculation of p. Assuming small coupling g/uo 1 of the 
first qubit and the resonator, we have then considered the TFIC as a linear bath for 
the resonator, and this approximation allowed us to calculate the resonator spectrum 
S in the coupled system. In this section, we add some remarks on the interpretation of 
the autocorrelator, the transition — > oo, and the linear approximation. Moreover, we 
discuss how a small amount of disorder in the qubit parameters due to imperfections in 
the fabrication process affects the resonator spectrum S. 
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Figure 2. Imaginary part of the qubit autocorrelator p(t) = (<7*(t)<7*(0)) of the 
transverse- field Ising chain with normalized transverse field £ = ft/ 2 J = 8 in the cases 
N = 20 (black) and N —> oo (magenta). 



3.1. Time- dependent correlations in the transverse- field Ising chain 



By means of the spin-free-fermion mapping described in section |2.3[ one readily finds 
P(t) = (^(t)al(0)} = J2^ itAk - (12) 

k 

Here and in the following, expectation values are calculated under the experimentally 
realistic assumption of zero temperature. In the uniform case Qj = Q and Jj = J ', 
where explicit expressions for 0^ and can be found, the limit iV — > oo can be taken 
analytically and yields [20J 

Pit) - 0(1 - |€|)(1 - + | l'M 1 + (,% k cosk e-«™. (13) 
Here, Q(x) is the Heaviside step function and A(k) stands for with continuous k 



[equation ( 10 )] . The first term on the RHS of ( 13 ) causes a nonzero mean value of Rep(t) 



in the ordered phase. Figure shows Imp(i) for £ = 8 in the cases = 20 [equation 



(12)] and A^ — > oo [equation ( [131 )] (the time evolutions of Rep and Imp are qualitatively 
similar and agree for |£| ^> 1 up to a phase). For small times, the curves coincide (the 
second covers the first). However, the finite size of the TFIC with A^ = 20 causes a 
revival of p at T r w 2N/v with v = max[dA(fc)/d/c] (v = 2J|£| for £ < 1 and v — 2J for 
|£| > 1). This can be understood in the following way. The autocorrelator p is related 
to the linear response A{o~l)(t) of the TFIC to a perturbation oc 5{t)a], relative to the 
equilibrium value (a]) = 0. Indeed, Kubo's formula predicts A(al)(t) oc Imp(t). The 
5-pulse at t — forces the first spin in the — x-direction. This local excitation in position 
space is composed of many excitations in A;-space. Since most of them have velocity v 
[20] . the local excitation propagates with velocity v through the system, is reflected at 
the far end of the chain, and causes revivals of p at multiples of T r = 2N/v. To further 
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clarify the transition N — > oo, we note that for large t, p has a standard deviation from 
its mean oc l/y/N. This can be expected from (12) since |p(t)| 2 ~ 1 / / iV s ''^2 k k , e 11 



and, for t — > oo, all terms in the sum except for those with k = k' will cancel. In general, 
the t — > oo fluctuations that we find for all time-dependent observables considered in 
this work are due to the finite system size and decrease with TV (but not all of them 
behave as oc l/y/N). 

3.2. Spectrum of the resonator - the linear approximation 



Taking the Fourier transform of equations (12) and (13) yields the spectrum p(u) of 
the TFIC for a force that couples to a\ for finite N and iV — > oo, respectively. In 
order to calculate the spectrum S of the resonator, whose coordinate (a) + a) couples 
to o\ [cf. equation (pfl)], we have suggested [20] a useful approximation: We consider 
the TFIC as a linear bath for the resonator. That is, we replace the TFIC by a set 
harmonic oscillators having the spectrum p of the TFIC. This approximation can be 
straightforwardly generalized to other contexts, where a different many-body system 
couples to a resonator. It is justified in the limit of small qubit-resonator coupling 
g/wo < 1, as we discuss in the following. 

The linear approximation for the TFIC-bath fails as soon as probing the resonator 
sufficiently excites the TFIC so that its nonlinearity becomes important. Thus, the 
linear approximation requires small coupling g and is worst if the TFIC is on resonance 
with the resonator (uo within the band of the TFIC). The "most nonlinear" bath 
possible for the resonator, that is, the bath whose nonlinearity becomes important for 
the smallest value of g, is a bath consisting of only a single qubit on resonance with 
the resonator. If the linear approximation is adequate for such a system in the limit 
g/ojQ <C 1, it will be also sufficient for our purposes. Therefore, we now consider the 
case N — 1 and Q = ujq of equation (IT]) and calculate the spectrum of the resonator by 
linearizing the single-qubit-bath. Since the atomic Hilbert space is small for N — 1, we 
can then numerically check the accuracy of our approximation. We also compare our 
approximation with the resonator spectrum calculated analytically within the rotating 
wave approximation (RWA), which is the standard approximation of % in this specific 
situation. 

For N — 1 and VL = cuq, the Hamiltonian H [equation (IT])] becomes 

Wn=i = uoa)a + g(o) + a)a x + L ^-a z (14) 

That is, the resonator coordinate (a^+a) couples to a single-qubit-bath with Hamiltonian 
1-L q = ^f-c z via a force ga x . The spectrum of this force is 

F q (co) = J dte^ q {ga x (t) ga x {0)) q = 2ng 2 5(cu - w ), (15) 
where the time evolution of a x and the expectation value q ( . ) q are to be calculated with 
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respect to (the ground state of) 7-L q . Now we linearize the system and replace Hn=i by 
Hun = u Q <Ja + g\a ] + a){b ] + b)+ wb ] b (16) 



with bosonic b,tf and parameters g' and w to be determined. In (16), the resonator 
couples to a force g'(b'+b) exerted by a bath that consists of a single harmonic oscillator 
with Hamiltonian l-i\ lo = wb^b. The spectrum of this force reads 

F ho (u) = 27i(g') 2 5(Lj-w). (17) 

Thus, we choose g' — g and w = u such that F ho 
calculate the autocorrelator of the resonator coordinate 

Piin(t) = im(W(t) + a(t)][a f (0) + a(0)]) Un , 
and its Fourier transform, the resonator spectrum 



F q . With this substitution, we now 



(18) 



Plin(Wj 



dte^Viin(t), 



(19) 



according to Equation (16). To that end, we express the resonator coordinate (a* + a) 
in terms of the (bosonic) eigenmodes c± with frequencies u± = a/ w o ± 2guo of T-Lvax, 



(a t + a) 




Using ( 20 ) , one readily finds 

Plin(t) 



2 



+ 



-iu)—t 



00- 



Plin(w) = 7T 



■^o(u; — w + ) + -^o(u; — u;_j 



(20) 



(21) 



(22) 



Before we go on and compare these approximate analytical results with numerical 
finite-size calculations for 1-Ln=i [Equation (14)], we calculate the same quantities on the 
basis of the standard approximation to T-Ln=i for g/<^o *C 1, the RWA (see, e.g., [32"]). 



This will be a helpful benchmark for estimating the quality of the linear approximation. 
In the RWA, the Hamiltonian T-Ln=i reduces to the Jaynes-Cummings Hamiltonian 

"Hrwa = UJ a^a + g(a^a~ + aa + ) + -^-a z . (23) 



This Hamiltonian can be straightforwardly diagonalized, and one can therefore 
analytically calculate the autocorrelator Prwa(^) and the spectrum prwa(^) of the 
resonator in the approximation provided by Hrwa, 

PrwaW = \ [e-^+a) + e""^"^] , (24) 
PhwaM = n{5(u - (u + g)) + 5(u - (uj - g))). (25) 



On the basis of (23), the results (24) and (25) are exact. 

The autocorrelator and the spectrum of the resonator can also be calculated 
numerically after truncating the photonic Hilbert space. This is achieved by expanding 
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Figure 3. Comparison of the rotating-wave approximation and the linear 
approximation with highly accurate finite-size numerics for a resonator with frequency 
wo resonantly coupled to a single qubit with coupling strength g/ojQ = 0.12. (a) 
Autocorrelator p(t) = ([a^(t) + a(t)][aJ(0) + a(0)]) of the resonator (red), and the 
same quantity calculated within the rotating-wave approximation (prwa, blue) and the 
linear approximation (pu n , green), (b) Spectrum p{u>) = J dte lujt p(t) of the resonator 
(red), and the same quantity calculated within the rotating-wave approximation (prwa, 
blue) and the linear approximation (prm, green). The dashed line is a guide to the eye. 



Hn=i an d the resonator coordinate (a^ + a) in the product basis {\s z , v)}, where s z =t, 4- 
and v G INo, and dropping all matrix elements with v > z/ max - In this finite-size 
approximation, the eigenvalues E n and eigenvectors \n) of "Hv=i can be numerically 
calculated (n = 0, . . . , n max = 2z/ max + 1) and give pit) and p{u) according to 

p(t) = J2 e- i(E "- Eo)t | (0| (a f + a) \n)\ 2 (26) 

n=0 

'Trnax 

~p{u) = 2nJ2 ~ (E n ~ E )) |(0|(at + a )\n)\ 2 (27) 

n=0 

Even for relatively strong coupling g/uo = 0.3, the numerical results for p(t) and p(u) 
are already converged if i/ max = 3 photonic excitations are taken into account. However, 
to be on the safe side, we choose z/ max = 10 in our calculations, which is still numerically 
easily tractable. 

Our results for the autocorrelator p(t) and the spectrum p(cu) of the resonator in 
%n=i are plotted, respectively, in figure |3]( a) [eqations (|2~T|),(24),(26)] and figure |3](b) 



[equations (22), (25), (27)]. In both plots, we choose g/u>o = 0.12, which is the largest 
ratio of g/ooo used in this work and in (2D]- The autocorrelator p(t) of the resonator 
(red) is well approximated both by the RWA (prwa, blue) and the linear approximation 
(pii n , green), and the quality of these approximations is essentially equal. For small 
t, the linear approximation might be even more accurate than the RWA, but becomes 
worse at large t. This can be understood in the frequency domain. In figure [3 
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we plot the spectral weights of the delta-peaks in the spectra p, Prwa, and pu n (red, 
blue, green) at the corresponding peak positions. The spectrum p contains also delta- 
peaks at higher frequencies than the ones plotted, but their weight is virtually zero 
(27r|(0|(at + a)|n)| 2 < 1CT 7 for all n 7^ 1,2). Both approximations yield good predictions 
for the positions and the spectral weights of the peaks in p. The RWA is more precise 
in predicting the peak positions and the linear approximation in predicting the spectral 
weights (note, however, that the peak positions in p lin and Prwa agree up to first order 
in g/oJo). Thus, the linear approximation is more precise for small t, in particular at 
t m and where the envelope of p(t) has a minimum, but becomes worse for large 
t. In summary, we conclude that even for the situation N = 1 and Q = uq, the 
linear approximation yields good results for the autocorrelator and the spectrum of 
the resonator in the limit g/ooo <C 1 that are qualitatively comparable to the usual 
RWA in this context. This implies that the linear approximation is well-justified in our 
calculation of the spectrum of a resonator coupled to a TFIC 



3. 3. Spectrum of the resonator - disorder effects 

The linear approximation for the TFIC allows one to express the spectrum S(u) of the 
(coupled) resonator as a function of the spectrum p(u) of the TFIC 



q(l A = 4e(u)[K + g 2 p(u)} 

1 ' [ w 2 M-wo-W)] 2 + [ k + 5 2 p»] 2 ' 1 ' 

Here, k is the full linewidth at half maximum of the Lorentzian spectrum of the 
uncoupled (g = 0) resonator and xi^ 2 ) denotes the principal- value integral xi^ 2 ) — 
l/(27r) J dQp(Q)Q/ (u 2 — Q 2 ). This result is actually general and holds for any linear bath 
coupled to a resonator, with an arbitrary spectrum p. Plots of S, with p{uj) being the 



Fourier transform of (13), are presented in [20]. However, in an actual implementation 
of the proposed setup, the qubit parameters Jj and Qj will not be perfectly uniform, 
due to imperfections in the fabrication process. We now investigate how this modifies 
the characteristic features of the spectrum S of the uniform system predicted in [20J. 
It is known in the field of random-matrix theory that disorder would have to be very 
strong in order to have a dominant effect on (average) spectra. We will observe the 
same here, in this concrete model system. 

For a nonuniform TFIC, no closed analytical expressions for p(u) are available. 
Thus, we have to consider finite system sizes and calculate numerically the relevant 
quantities, specifically, the spectrum of a finite-size nonuniform TFIC, 

p(u) = 2ir^l 1 6(u-A k ), (29) 



which is the Fourier transform of Equation (12). To take the effect of qubit decay 



processes into account, we phenomenologically broaden the delta-peaks in (29) and 



replace them by Lorentzians of width 7 around the A*. We model the nonuniformity 
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Figure 4. (a) Spectrum S of a resonator coupled to a finite uniform TFIC with 
N = 20 vs. probe frequency w and normalized transverse field £ = Q./2J. The 
parameters are g = 0.12, J = 0.08, k = 10~ 4 , and 7 = 5 x 10~ 3 (in units of lu )- 
For better visibility of the features, values > 3 are plotted white, (b) Spectrum S(ui) 
for £ = 6.1. This curve corresponds to a cut along the arrows in (a), (c) Same as in 
(a) but with Qj and Jj following a Gaussian distribution with a standard deviation of 
2% around their mean values, (d) Cut along the arrows in (c). 



of the qubit parameters by writing Qj = Qtj and Jj = Jr'j and choosing Tj and rj to 
be random variables, which follow Gaussian distributions with means 1 and standard 
deviations o T = o T i = 0.02. With a typical set of system parameters that was also used 
in [20], we numerically calculate p(u) according to (29) and the corresponding resonator 
spectrum 5* according to (28 ). In order to judge the effects of disorder, we also reproduce 



our calculation of S for the corresponding uniform system [20] (Fig- S6). Figure [4] shows 
S as a function of u and the (mean) normalized transverse field £ = Q/2J for the 
uniform system [figures |i](a,b)] and for a typical disorder configuration [figures |4]^c,d)]. 
In the uniform case, the signatures of the QPT at £ = 1, the dispersive shift of the 
resonator frequency, and, on resonance, the double-peak with a separation of 4 J (rather 
than 2g as in the case N = 1) that we have discussed in detail for iV — > 00 in [20j are 
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clearly visible also for N = 20. These characteristic features are insensitive with respect 
to a small amount of disorder in the system parameters, as figures (4tc,d) demonstrate. 

4. Disorder effects on the system dynamics 

A particularly interesting application of the proposed system would be to simulate the 
nonequilibrium dynamics of the TFIC. In [20], we have suggested to experimentally 
track the propagation of a localized excitation in the (uniform) TFIC that can be easily 
created in our system and to measure the system dynamics after quenching the transition 
frequencies of all qubits. In this section, we show that none of the predicted experimental 
results changes qualitatively if the parameters of the TFIC are slightly disordered, as has 
to be expected in reality. However, strong disorder, accessible by deliberately detuning 
individual qubits, is shown to produce qualitatively different physics in the experiments 
proposed in [20] , like Anderson localization of the propagating excitation. We develop a 
quantitative theory of the disorder effects on the system's nonequilibrium dynamics that 
explains the results of numerical experiments with the disordered TFIC. This theory 
might be helpful for experimentalists to estimate system and disorder parameters for 
successfully performing nonequilibrium experiments with the TFIC without having to 
do numerical simulations. 

4-1. Propagation of localized excitations 

For the first type of experiments we have suggested in [20] , it is assumed that the TFIC 
is deeply in the paramagnetic phase (£ ^> 1) and detuned from the resonator. In this 
situation, the TFIC is essentially decoupled from the resonator and its ground state is 
characterized by (a 3 z ) ~ —1. Applying a fast 7r-pulse to the first qubit thus creates a 
localized excitation in the system that subsequently propagates trough the chain due to 
the qubit-qubit coupling J . The time evolution of the observable (a 3 x ) after the 7r-pulse 
can be approximately described by [20] 

W>(*) = - I>*A* + 5> i(Afc - Afe,)t [0 fc ,i^',i 

b b b 1 

x (ipk,j4>k',j + ipk>,j(f>k,j)]- (30) 

We plot this result in figure |5](a) for all j in a chain of length N = 20 and for a 
mean normalized transverse field £ = VL/2J = 8 (same system parameters as in [2D]), 
and we randomly choose Qj = Qtj and Jj = Jt'a according to Gaussian distributions 
with standard deviations of 2% from the mean values Q and J as before (right panel). 
The experimentally measurable observable (c*)(t) is singled out in the left panel. The 
propagation of a localized excitation through the chain, and its reflection at the far end 
of the chain that leads to a distinct revival of (crl)(t) at t ~ Nj J [20], are still clearly 
visible in this slightly nonuniform system. 




Figure 5. Propagation of a localized excitation in a weakly disordered transverse- field 
Ising chain of length N = 20. Specifically, the density plots show the nonequilibrium 
time evolution of (a 3 z ) for all j after a 7r-pulse on the first qubit while the system is in the 
paramagnetic phase (mean normalized transverse field £ = 8). For better visibility of 
the features, values > —0.5 are plotted white. The experimentally accessible observable 
(al) is singled out in the left panels, (a) The qubit transition frequencies flj and 
qubit-qubit couplings Jj are randomly chosen according to Gaussian distributions 
with standard deviations of 2% from the mean values, (b) Same as in (a) but with 
qubit f 1 strongly detuned. 



If the transition frequencies Qj of the qubits can be tuned individually, the effective 
length of the TFIC has been shown to be adjustable by strongly detuning one qubit 
from the others [20]. This holds true also for a slightly nonuniform system: Figure |5|b) 
shows a typical result for a system with the same parameters and disorder strength as in 
(a), but with qubit 11 strongly detuned by setting m = 1-3. This result is qualitatively 
identical with the result for the corresponding non- disordered system J2U]. The strong 
nonuniformity at j = 11 acts as a barrier for the propagating excitation and leads to its 
reflection. Thus, it effectively changes the length of the TFIC. 

Having shown that the experiments with propagating localized excitations proposed 
in [20] yield qualitatively the same results for ordered and slightly disordered systems, 
we now proceed and study disorder effects on this type of experiments quantitatively. 
Parts of the following analysis also apply to other nonequilibrium experiments with the 
TFIC, as will be discussed in the context of the quantum quenches (section 4.2). 

Since it is assumed that the system is deeply in the paramagnetic phase, the 
mean qubit transition frequency Q is larger than the modulus of the mean qubit-qubit 
coupling J, Q/J 1. As before, we further assume uncorrelated disorder of the system 
parameters via Qj = Qtj and Jj = JtL where Tj and rj follow Gaussian distributions 
with standard deviations o T and a T i from 1. That is, for a T = cr T i, the absolute variation 
of the Qj will be larger than the absolute variation of the Jy Therefore, the dynamics 
of the system may be expected to be much more sensitive to increasing a T than a T < . 
Moreover, one may expect that disorder effects start to qualitatively affect the system 
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dynamics even of small systems (that is, on the scale of neighbouring qubits j and 
j + 1) when the disorder in the qubit transition frequencies becomes comparable to the 
modulus of the mean qubit-qubit coupling, Vta T = J. These deliberations are confirmed 
by numerical experiments: We first consider the wave functions gkj and hk,j in position 
space (rjk = Ylj=i9k,j c j + ^fcj c ])- F° r zero disorder, they are extended over the whole 
chain (except for the mode with — > in the ordered phase |33]). Increasing a T 
localizes the wave functions much more strongly than increasing a T >, and the localization 
length of the wave functions indeed reduces from many (3> 1) sites to a few (> 1) sites 
at Qa T ~ J. Correspondingly, the propagation of an excitation initially localized at site 
1 is only weakly affected by disorder in J . However, if a T > J/Q, it propagates only a 
few sites before becoming completely trapped due to the disorder. This manifestation 
of Anderson localization [31] is illustrated in figure [6|a), where we have used the same 
system parameters as in figure [5j but we have randomly chosen Tj and rj according to 
Gaussian distributions around 1 with standard deviations cr T = J/Q = 0.0625 = a r i. 
For definiteness, we always choose o T > = a T in the following. 

We have seen that for |£| ^> 1 (paramagnetic phase) the effective disorder strength 
in the quantum Ising chain is set by a T Q/J oc ov|£|. Now we try to determine how 
the relevant observables in the currently considered type of experiments depend on 
this quantity. The observable we focus on in the following is the maximum excitation 
probability (maximized over time) of the jth qubit caused by the propagation of the 
localized excitation through the disordered chain. In an experiment, one would for 
instance create an excitation of the first qubit and measure the excitation probability 
of some other (e.g. the A^tli) qubit as a function of time. The maximum excitation 
probability of the jth qubit is an important quantity since it will determine if the effect 
of the propagating excitation can be measured at site j, given a certain measurement 
resolution. In a single disordered system, the maximum excitation probability of qubit 
j will depend on the specific (random) disorder configuration of this system. Therefore, 
a study of the effect of disorder as characterized by the statistical quantity u T \C\ can 
only refer to the statistical average of the maximum excitation probability of qubit 
j in one disordered system over an ensemble of many disordered systems (disorder 
configurations), all chosen according to the same probability distribution. Stated as a 
formula, this ensemble average of the maximum excitation probability of qubit j is given 
by 

pi T1 i f i = ^(mp[<4M+i). (si) 

Here, the double overbar ~ denotes the ensemble average over many disordered systems 
(disorder configurations) with the same system and disorder parameters £, J, and 
a T = o T i . This average is taken after one has maximized (cr^)(£) for a specific disordered 
system over time. Our goal is to find the explicit functional dependence of p 3 ^ on 
a T and |£| (in fact, we expect dependence only on the product oy|£|). Note that we 
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assume that p J a ^ depends neither on the sign of £ nor explicitly on the mean qubit- 
qubit coupling J, but only on the ratio of Q and J (via |£|). This is strictly true for 
a T = o T i = 0. By explicitly solving equations ([6| for this case (20], one can show that 
after substituting £ — > — £, the new allowed wave vectors are q = 7r — k with A g = A*, 
= {—l) N ~i(j)kj, and = (— ^) N ~^k,j- With that one can easily see that equation 



(30) does not depend on the sign of £. Moreover, 0^ and ipk are independent of J [which 
also follows from equations (]6])], and A& oc J such that changing J corresponds only to 
a rescaling of time. The influence of disorder, however, is essentially set by cr r |£| (for 
|£| 3> 1), as we have argued above. Consequently, we may take p> to be independent 
of J and of the sign of £. Nevertheless, to keep notation short, we write £ instead of 
|£| for the remainder of this chapter. For simplicity, we first focus on a semi- infinite 
system (N — >■ oo) and discuss later the increase of p- 7 at the end of the chain (due to 
the refocusing of the dispersed wave packet of the propagating excitation). 

As usual for disordered systems (e.g. [35]), we will try to characterize the disorder 
effects on the ensemble-averaged maximum qubit excitation p 1 ^ g via a mean free path. 
To that end, it pays to first discuss in more detail the uniform Even 
there, analyzing the propagation of the dispersive wave packet that determines the 
maximum excitation probability of a qubit requires some care. For £ ^> 1, this 
excitation probability does not depend on £. This is because the dispersion relation 
of the TFIC becomes that of the tight binding model, A^ = 2Jyl + £ 2 — 2£ cos k ~ 
2Jsign(£)(£ — cos k). Thus, £ only sets the band gap but does not influence the shape of 
the dispersion relation. Except for the aforementioned boundary effects, p 3 Q ^ also does 
not depend on N. This is evident from Figure 6^c), where we plot pj^ for several £ and 
N. The curves for different £ but same N lie almost on top of each other (henceforth, 
we drop the index £ form p^f), and curves for different can be distinguished only 
by the boundary effects, that is, by the strong increase of p J at j = N (which will be 
discussed later). The decay of the p J with j is relatively slow (slower than 1/j), which 
should considerably simplify the experiments proposed in [20]. This slow decay of p J Q 
can be understood from the dispersion relation A& of the system which, for £ 3> 1, is 
quadratic in k at k ~ 0, tt, and linear at k ~ 7r/2: If an initially localized wave packet 
with width s and momentum q, ip(x,0) = ae~ x2 ^ 2s2+tqx , a = (s 2 7r) -1 / 4 , is evolved in 
time by the Hamiltonians Hi = h\h and H2 = h,2k 2 , respectively, one finds 

mx,t)\ 2 Hi = |^(x-M,0)| 2 = aV^ 2 / s2 , (32) 
2 a 2 s 2 ( (x — 2h 2 tq) 2 



^'"^^•"l-wJ- (33) 

That is, for Hi, maximum and width of the probability distribution for finding the 
particle at a position x are constant, while for H 2 and strong initial localization (or 
large times) the width is cx t and the maximum ex 1/t. As the dispersion relation of the 
TFIC interpolates between these two cases, one may expect a decay of p^ slower than 

Vi- 



CONTENTS 



19 




Figure 6. (a) Propagation of an initially localized excitation in a strongly disordered 
transverse-field Ising chain. Initialization and system parameters are identical to 
figure [5]j a), but the Qj and Jj are randomly chosen according to Gaussian distributions 
with standard deviations of 6.5% from the mean values. The plot clearly shows that 
strong localization of the excitation prohibits its propagation through the chain, (b) 
Mean free path I of the propagating excitation (defined in the main text) vs. normalized 
standard deviation a T of the qubit transition frequencies for different values of the 
normalized transverse field £ on log-log scale. The points are c as gained by 
numerically averaging many disorder configurations. The lines are best fits of l/cr°£ b 
to these data. (c,d) Behavior of a non-disordered system, with uniform flj = ft and 
Jj = J , for comparison, (c) Maximum excitation probabilities p 3 * of the jth qubits 
in the nonequilibrium time evolution of uniform transverse-field Ising chains of lengths 
N = 10, 20, 30, 40, 50 after the first qubit has been flipped. For each chain length, p J ^ 
is plotted for £ = 3,5,8 (red, green, blue). Apart from boundary effects, the decay 
of j»g £ with j is slower than oc 1/j. The maximum excitation probabilities p^^ of the 
last qubits of the chains are significantly enhanced compared to nearby bulk sites, (d) 
Maximum excitation probability p(f of the ATh qubit vs. chain length N (for any 
£ » 1). 



Coming back now to the disordered case, one might suspect that the ensemble- 
averaged maximum qubit excitation p^ g is related to the corresponding quantity for a 
non-disordered system p J via an exponential decay, governed by a finite mean free path 
l aTt £ for the propagation of the localized excitation, 

pL T t=P' e- j/l '>T>t. (34) 



If (34) holds, 




pJ ° ; (35) 

should be independent of j. This observation can be used to check our assumption 



(34). We numerically calculate g for all combinations of £ = 3, . . . , 8 and 100 x u T £ 
{1, 1.5, 2, 2.5, 3, 4, 5, 6, 7, 8} in a chain of length iV = 20, and we average over 100 disorder 
configurations. This turns out to be a good compromise between calculation time and 
ensemble and system size as long as the effective disorder cr T £ is not too small (see 
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below). With these p^ we calculated the RHS of (35) for j — 5, . . . , 16. Other j are 
not considered, in order to minimize boundary effects. Our results for j = 5, . . . , 16 are 
approximately equal, with the ratio of standard deviation to mean value being < 0.1 
for given a T and £. We note that for very weak effective disorder cr T £ < 0.1 we have 
to average over 500 disorder configurations such that this ratio is < 0.1, because with 
decreasing ratio Po/p^ g the slope of the logarithm on the RHS of (35) increases. Thus, 
the numerical data seem to confirm our assumption (34), and the influence of disorder 
on the considered experiment is captured by a mean free path l aT £. In our subsequent 
analysis, we try to find simple expressions for this quantity. 

The propagation of the localized excitation in the Gaussian disordered TFIC is 
akin to the propagation of a particle in an uncorrelated random potential V(r) with 
(V(r)V(r')) = Vq5(t — r'). To lowest order in perturbation theory (Fermi's golden rule, 
e.g. [35]), the mean free path of the latter decreases as the inverse square of the disorder 
strength, oc 1/Vq 2 . In our case, the effective disorder strength is determined by the 
dimensionless quantity <7 r £. Therefore, we expect that 

t*rA = j^yr (36) 

To check this, we calculate l aT g for the same combinations of o T and £ as before by 
averaging the RHS of (35) over j = 5, . . . , 16. Then we fit the function l(cr T , £) = l/u^ b 
to our data for 1 Ut £- We find the exponents a ~ 2.002 and b ~ 2.071, which comes 
close to our expectation of a = b = 2. Numerical data and fit are plotted on log-log 
scale in figure |6](b). As long as the effective disorder strength is not too big (<7 T £ < 0.2), 
l(cr T ,£) with the fit values of a and b reproduces the numerically (by ensemble-averaging) 
extracted mean free path l ar £- Here, one may attribute the deviations of a and b from 2 
to the finite ensemble sizes. For stronger disorder, however, the fit of l(cr T ,£) = l/<r"£ fe 
begins to deviate from l„ T ^. Thus, higher-order effects (beyond Fermi's golden rule) 
and/or the disorder in J seem to be no longer negligible. 

Finally, in setups with a second readout resonator (cf. figure [I]), the maximum 
population p^ T ^ of the iVth qubit will be an experimentally relevant quantity. Since 
the dispersed wave packet of the propagating excitation is refocussed at the end of the 
chain, the maximum excitation probability of the iVth qubit is considerably enhanced 
compared to nearby bulk qubits [see figure |6](c)]. It turns out that p^ T £ can also be 
estimated by means of (34), the mean free path (36), and the value of p$ for the 
corresponding non- disordered system, which we plot for iV = 1, . . . , 50 in Figure |6](d). 

Summing up, equations (34) and (36), together with Figures [6]^c,d) allow one to 
easily estimate suitable system and disorder parameters for successfully implementing 
the presently considered type of experiment. For instance, if in a system with iV = 30 
and £ = 4 the iVth qubit should get a population of p^ T ^ = 0.3 (which corresponds 
to max [(cf)(t)] = —0.4), then one can roughly (i.e., averaged over many systems) 
afford a standard deviation of the qubit transition frequencies from their mean of a T = 
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qubit number j (a z ) qubit number j (°z ) 



Figure 7. Propagation of a localized excitation in a nonuniform transverse-field 
Ising chain of length N — 30 and with normalized transverse field £ = 4. (a) The qubit 
transition frequencies flj and qubit-qubit couplings Jj are randomly chosen according 
to Gaussian distributions with standard deviations of 3.4% from the mean values, (b) 
Same as in (a) but with a cosine modulation of the qubit transition frequencies flj 
with standard deviation w 3.4%, instead of uncorrelated disorder of flj and Jj. 



[(iV^ 2 )" 1 log (pf/0.3)} 1 ^ 2 0.034, where we have extracted pl° 0.53 from figure 6^d). 
A typical result for these parameters is plotted in figure fna). Here, the maximum 
excitation probability is found to be Po°034 4 ~ 0-29 (since max t [(crf°)(t)] ps —0.43). 

We remark that the foregoing deliberations only hold for uncorrelated disorder of 
the system parameters and do not take into account qubit decay. Correlated disorder 
can yield qualitatively different results and has to be studied explicitly via equation 
(30). We also remark that if the Qj are individually tunable, it becomes possible to 
study the propagation of localized excitations in arbitrary potentials. For instance, it 
might be interesting to choose Qj = + y/2a T cos(2nj /N)] and to compare the system 
dynamics with the Gaussian disordered case. For large N, both distributions of Qj have 
the same mean and the same standard deviation, but in the former case the system is 
not disordered and the localization of the propagating excitation is much weaker than 
in the genuinely disordered case. Figure [7|b) shows the propagating excitation in such 
a system with N, £, and a T as in figure IWa) (with uniform Jj). 



4-2. Quench dynamics 

The second type of nonequilibrium experiments we have proposed in [20] relies on the 
possibility to rapidly change the transition frequency Q of a superconducting qubit in 
a circuit QED system by tuning the magnetic flux through its SQUID loop. This has 
been shown to be possible virtually instantaneously on the dynamical time scale of a 
circuit QED system (H El [7] , without changing the system's wave function. Let us now 
assume that the circuit QED quantum simulator of the (uniform) TFIC proposed in [20] 
is implemented with Cooper-pair boxes. For this system, such a sudden sudden change 
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of all Vtj = Vt corresponds to a global quantum quench of the normalized transverse 
magnetic field £ = VL/2J. We remark that one can also produce quenches of £ by 
using transmons in a non-standard parameter regime instead of Cooper-pair boxes, or 
by using usual transmons with tunable coupling capacitances (36J, [20] . 

The nonequilibrium dynamics of the uniform TFIC following a quantum quench is 
currently subject to much theoretical research, e.g. [221 EH EE1 ESI HQl SI] , and should be 
experimentally observable with our circuit QED quantum simulator. In this context it 
is usually assumed that for t < the system is in the ground state |0) a of a Hamiltonian 
Hi a (characterized by £ a ). At t — 0, the overall transverse field is changed, £ a — > 
and the nonequilibrium time evolution of some observable O under Hij, is investigated, 

(0)(t) = a (0\e m ^Oe- m ^\0) a . (37) 

In [20] we have focussed on the time evolution of the local transverse magnetization 
{a{) and the end-to-end correlator {cr].a^) (indicating long-range order) after quenching 
£ within the paramagnetic phase. These quantities should be experimentally easily 
accessible in our system. In this section, we show that also for such quantum quenches 
the predicted experimental results of [20] are insensitive to a small amount of fabrication- 
induced disorder. 

In general, two sets of the Q- and ^-parameters, {O" } and {jf h }-, fully specify 
the Hamiltonians Hi A /b [equation Q2p] . Given these parameters, the time evolution (37) 
of the local magnetization and the end-to-end correlator can be written as [20J 



b k',j 

k,h' 

b i A 6 \ i \r „„„ jV A b 



[X Kkl cost{k b k + A b fc ,) + y M ,cost(A* - A*,)]}, (38) 
[X fcjfc/ cost(A^ + At) - n )fc ,cost(A£ - A^]}. (39) 



Here, 



X k , kl = [{g b k ) T H a + (h b k ) T G a ] \{G a ) T g\i + (H a fh b k ,], (40) 
Yh» = i(9 b k ) T H a + (h b k ) T G a ][(H a ) T g b f + (G a ) T hU, (41) 

and G and H are matrices that respectively contain the g k and h k as columns. In these 
equations, a quantity carrying the index a or 6 is to be calculated from equations ^ 
with parameter set a or b. 

To implement disorder of the system parameters before the quantum quench, we 
write again = fi°Tj and Jj = Jt'^ and we randomly choose Tj and rj according to 
Gaussian distributions with standard deviations a T and oy from 1. As we have argued 
in section |2~7L tuning the flux $ through the SQUID loops of the qubits only changes 



the mean qubit transition frequency Q a — » fi b (and, thus, the mean transverse field 
£ a = VL a /2J —>•£& = VL b /2J), but leaves r,-, J7, and r- unaffected. Hence, by fixing £ a / fe 
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Figure 8. (a) Time evolution of the magnetization (cr^) in a disordered TFIC of length 
iV = 30 after a quench of the mean normalized transverse field £ = £1/2 J = 8 — >• 1.2 
(right). Values < —0.9 (> —0.6) are plotted black (white). The measurable observable 
(al) is plotted separately in the left panel (black), along with the corresponding trace 
for a uniform system (green), (b) Time evolution of the end-to-end correlator (c^a^) in 
a disordered TFIC of length N = 30 after a quench of the mean normalized transverse 
field £ = 8 — > 1.5 (black), along with the corresponding trace for a uniform system 
(green). In both plots the qubit transition frequencies flj and qubit-qubit couplings 
Jj are randomly chosen according to Gaussian distributions with standard deviations 
of 2% from the mean values f2 and J . 



and a T / T ', the system is fully specified before and after the quench (as in section 4.1 , the 
absolute values of fW 6 and J can be absorbed in the time scale Jt of the dynamics), 
and we are ready to evaluate equations (38) and (39). 

Figure pB[a) shows the local magnetization (cr 3 z )(t) for all j and for the same system 
parameters as in figure 4 of |20j, but with fL and Jj having standard deviations 
a T = o T i = 2% around their mean values (right panel). The experimentally easily 
measurable trace of (al) is singled out in the left panel (black). For comparison, we also 
plot (green) the local magnetization of the first qubit (al) of the uniform system (as 
plotted in the left panel of figure 4 of [2D])- Correspondingly, Figure |8](b) shows Equation 
(39) for a uniform system as in Figure S6 of [20] (green), and with 2% disorder in Q^ b 
and Jj (black). The plots demonstrate that the quench dynamics of the considered 
observables is not qualitatively affected by the presence of a small amount of disorder. 

For a more systematic analysis of the disorder effects on the quench experiments 
considered here we make use of our findings for the mean free path of a propagating 
localized excitation from the previous section. This is possible because the quench 
dynamics of the TFIC is governed by the propagation of quasiparticles (QPs) through 
the system [38] [39j SOI HU H2, [20] . These correspond to flipped spins, essentially like 
the localized excitation of the previous section. Indeed, if the system is initially in 
the paramagnetic phase, the time evolution immediately after the quantum quench 
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e -it?4|Q^ a j-j , e -itJ(€b-ta,)/£a.cri<4: + |q^ gjpg p a j rs Q f adjacent spins so that they point 
in the +z-direction. Due to the qubit-qubit coupling, these local excitations propagate 
as QPs with velocity v ~ 2J through the chain. For an interpretation of the quench 
dynamics and the time scales indicated in the plots in terms of these QPs, see [20]. If 
£b is in the paramagnetic phase, the mean free path I of the QPs in a disordered TFIC 
can be estimated by / = 1/ (er r £b) 2 according to the previous section. The characteristic 
quasi-T-periodic behavior (T = N/v) of the local magnetization after the quench in 
the non-disordered TFIC can be understood as a revival of coherence each time QPs 
initially generated at the same spot meet again [4T| [20] . This happens when the QPs 
have travelled multiples of the chain length N. If there should be a significant probability 
that two contiguously generated QPs meet again at least once before being scattered 
and thus decrease the local magnetization at t — T, the mean free path has to be 
sufficiently large, I > 2N. The appearance of significant end-to-end correlations after 
the quench (that are stronger than those for t — > oo) requires that QPs generated in the 
middle of the chain reach the edges of the chain without being scattered, hence I > N. 
We have performed numerical experiments which indeed suggest that the corresponding 
values of a T mark the transition from weak to strong disorder, where the described 
phenomena are no longer present. In that sense, the distinctive features of the quench 
dynamics of the end-to-end correlator are less sensitive to disorder than those of the 
local magnetization (and, due to the shorter time-scale, less sensitive to qubit decay). 
We finally note that also here the effective chain length can be adjusted by strongly 
detuning individual qubits (this can also be used to create local quantum quenches by 
'joining' two initially independent chains) and arbitrary effective potentials Qj can be 
chosen. 

5. Conclusion 

In the quest for controllable large-scale quantum systems, the framework of circuit QED 
offers several advantages, such as fast, high-fidelity readout, great flexibility in design, 
and steadily increasing coherence times. However, a potential significant disadvantage 
arises from the hardly avoidable static noise and disorder sources in these man-made 
devices. The central result of the present work is that also in this respect, there is reason 
to be optimistic: The requirements on the homogeneity of the system parameters for 
observing interesting (and predictable) many-body physics in a circuit QED system 
are not too high to be achievable with present-day or near-future technology. This 
underlines the prospects of circuit QED as a promising platform for implementing 
quantum simulations of complex quantum many-body Hamiltonians. In addition, we 
have shown that circuit QED quantum simulators could be used to study deliberately 
the effects of tunable disorder on quantum many-body dynamics. 
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